
% Function to calculate maximum rate of climb
%
% h = Altitude [m]
% W = Weight [N]
% Throttle [%]
% flag = 1 for jet
% flag = 0 for turboprop
% n = number of engines
%
% Tav        = Available Thrust [N]
% RC_max     = Maximum rate of climb [m/s]
% tas_RC_max = True air speed of maximum rate of climb [m/s]

function [RC_max, tas_RC_max] = calc_max_RC(h, W, Throttle, ISA, flag, n)

global ac oswald

S = ac.geo.wing.S;
AR = ac.geo.wing.ar;
eps = 1e-2;

[p,rho,mu,T,sound,delta_isa] = isa_cos(h,ISA,1);

% Jet
if flag == 1
    tas = 0:10:0.8*sound;
    
    RC_max = 0;
    tas_RC_max = 0;
    
    for i = 1:length(tas)
        Cl = W/(0.5*rho*tas(i)^2*S);
        Cd0 = calc_Cd0_aircraft(ac,h,tas(i));
        Cd = Cd0 + Cl^2/(pi*AR*oswald);
        
        [TSFC,Tav(i),FC] = engine_model(ac,h,tas(i),Throttle,ISA);
        Tav(i) = Tav(i)*n; % Thrust for n engines
        
        RC(i) = tas(i)*(Tav(i)/W - Cd/Cl); % Roskam Performance pg.383 Eq.9.28
        
        if RC(i) > RC_max
            RC_max = RC(i);
            tas_RC_max = tas(i);
        end
    end
    
% Turboprop
elseif flag == 0
    
    % Loop to find tas_RC_max
    tas = 100;
    Cd0 = calc_Cd0_aircraft(ac,h,tas);
    Cl = sqrt(3*Cd0*pi*AR*oswald); % Cl for best rate of climb. Roskam performance pg.399
    tas_aux = sqrt(2*W/(rho*S*Cl));
    while abs(tas-tas_aux) > eps
        tas = tas_aux;
        Cd0 = calc_Cd0_aircraft(ac,h,tas);
        Cl = sqrt(3*Cd0*pi*AR*oswald);
        tas_aux = sqrt(2*W/(rho*S*Cl));
    end
    
    % Roskam performance pg.399
    tas_RC_max = tas;
    Cd = Cd0 + Cl^2/(pi*AR*oswald);
    Drag = W*Cd/Cl; % [N]
    THPreqd = Drag*tas_RC_max; % Thrust HorsePower [W]
    
    [TSFC,Tav,FC] = engine_model(ac,h,tas_RC_max,Throttle,ISA);
    Tav = Tav*n; % Thrust for n engines
    THPav = Tav*tas_RC_max;
    
    RC_max = (THPav - THPreqd)/W; % [m/s]
end

